Identifying phase synchronization clusters in spatially extended dynamical systems 
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We investigate two recently proposed multivariate time series analysis techniques that aim at 
detecting phase synchronization clusters in spatially extended, nonstationary systems with regard to 
field applications. The starting point of both techniques is a matrix whose entries are the mean phase 
coherence values measured between pairs of time series. The first method is a mean field approach 
which allows to define the strength of participation of a subsystem in a single synchronization cluster. 
The second method is based on an eigenvalue decomposition from which a participation index is 
derived that characterizes the degree of involvement of a subsystem within multiple synchronization 
clusters. Simulating multiple clusters within a lattice of coupled Lorenz oscillators we explore 
the limitations and pitfalls of both methods and demonstrate (a) that the mean field approach is 
relatively robust even in configurations where the single cluster assumption is not entirely fulfilled, 
and (b) that the eigenvalue decomposition approach correctly identifies the simulated clusters even 
for low coupling strengths. Using the eigenvalue decomposition approach we studied spatiotemporal 
synchronization clusters in long-lasting multichannel EEG recordings from epilepsy patients and 
obtained results that fully confirm findings from well established neurophysiological examination 
techniques. Multivariate time series analysis methods such as synchronization cluster analysis that 
account for nonlinearities in the data are expected to provide complementary information which 
allows to gain deeper insights into the collective dynamics of spatially extended complex systems. 

PACS numbers: 87.19.La, 05.45.Tp, 05.45.Xt, 05.10.-a 



I. INTRODUCTION 



Spatially extended complex dynamical systems may be 
thought of being composed of numerous constituents (dy- 
namically formed subsystems) each having its own dy- 
namics. Typically the relevant state variables of such sys- 
tems can not be observed directly but only through some 
observation function that projects the high-dimensional 
state space onto an observation space of much lower di- 
mension, resulting in a set of time series. Multivari- 
ate analyses of such time series might then help to gain 
deeper insights into the collective dynamics of spatially 
extended systems. Although a number of time series 
analysis methods have been developed over the past (see 
[1-5] for an overview) most techniques either allow to 
characterize single time series (univariate approaches) or 
to investigate relationships between two time series (bi- 
variate approaches). However, applying bivariate tech- 
niques to pairs of time series - taken from a multichan- 
nel recording - does not necessarily allow to identify the 
relevant information in the full data set. The latter is of 
particular interest for scientific fields investigating spa- 
tially extended dynamical systems, such as meteorology, 
economics, social science or neurosciences, where a com- 
plex but relatively sparse connectivity between subsys- 
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terns prevails. Understanding brain function - both dur- 
ing physiological and pathophysiological conditions (as 
e.g. in the case of epilepsy) - requires a characterization 
and quantification of the collective behavior of neural net- 
works generating signals at different areas. 

In principle, multivariate time series analysis tech- 
niques can be used to investigate mutual relationships 
between arbitrary numbers of time series. A large va- 
riety of methods [6] aim at revealing additional infor- 
mation by classifying time series into different groups. 
In addition to the classical principal component analysis 
(also known as Karhunen-Loeve transform) [7] indepen- 
dent component analysis [8] provides a decomposition 
of data into independent source signals, and if the as- 
sumption of independence holds, it can be regarded as a 
suitable method. If independence can not be assumed, 
mutual information based methods might be more ap- 
propriate [9, 10]. The partial coherence [3] measures the 
fraction of coherence between two time series that is not 
shared with a third time series. Whereas the partial co- 
herence is based on the assumption of linearity and thus 
does not capture nonlinear interactions, the recently pro- 
posed concept of partial phase synchronization [11] was 
designed to account for nonlinearities of the dynamics 
under investigation. In order to study causal relations 
among simultaneously acquired time series generated by 
linear stochastic systems the Granger causality [12] can 
be used by fitting autoregressive models. Besides a re- 
cently suggested nonlinear extension of Granger's ideas 
[13], we mention the directed transfer function that is 
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defined for an arbitrary number of channels [14] and is 
based on a multivariate autoregressive model approach 
[15]. In Refs. [16, 17] the partial directed coherence has 
been proposed for inference of linear Granger causality 
in the frequency domain based on vector autoregressive 
models of appropriate order. Both methods have been re- 
peatedly applied to study interdependencies and causal 
relationships among neural signals (see e.g. [18-21] and 
references therein). 

Over the last decade time series analysis techniques 
known from random matrix theory [22, 23] have been 
repeatedly shown to allow an improved characterization 
of complex spatiotemporal correlation patterns. In these 
studies particularly the equal-time correlation matrix has 
been analyzed that was constructed from multivariate 
data sets obtained empirically in scientific fields rang- 
ing from economy [24, 25] and meteorology [26] to the 
neurosciences [27-31]. More recently an approach for the 
detection of clusters in financial data based on properties 
of the eigenvalue spectrum of the equal-time correlation 
matrix was proposed in Ref. [32]. By filtering out the 
random part and collective market-wide effects the au- 
thors were able to detect groups of stocks by optimizing 
the matrix representation. Such an approach, however, 
requires a clear-cut definition of the random part, which 
is usually assumed to be associated with the small bulk 
eigenvalues. In Ref. [29] it was demonstrated though 
that the lower part of the eigenvalue spectrum may con- 
tain essential information. 

Another way to study interactions in spatially ex- 
tended systems is based on a statistical analysis of phase 
synchronization phenomena [4, 33]. In most studies, 
however, the analysis of empirical multivariate data has 
been accomplished by a repeated application of bivari- 
ate synchronization measures and it remains to be es- 
tablished whether this approach allows to fully char- 
acterize a common integrating structure that may be 
present in the data. Addressing this issue, Allefeld and 
Kurths proposed a genuinely multivariate phase syn- 
chronization analysis method [34, 35] and successfully 
applied their method to electroencephalographic (EEG) 
data recorded during a psychological experiment. The 
authors concluded their method to provide additional in- 
formation on brain dynamics in a topographically, tem- 
porally, and frequency-specific way, as well as in other 
fields concerned with multivariate oscillatory processes. 
More recently, Allefeld and colleagues [36] introduced a 
new approach that addresses a limitation of their origi- 
nal method, namely the assumption of the existence of 
a single synchronization cluster in the data. By using 
methods from random matrix theory the resulting ap- 
proach appears to be capable of identifying multiple syn- 
chronization clusters in the data which makes it highly 
attractive to characterize pathophysiological, spatiotem- 
poral synchronization phenomena in multichannel EEG 
recordings from epilepsy patients. 

Recent findings indicate that bivariate analysis tech- 
niques allow to characterize physiological and pathophys- 



iological phenomena in the human brain (see Ref. [37, 38] 
for an overview) , and it can be expected that multivariate 
phase synchronization analysis techniques provide com- 
plementary information. To address this issue we here 
study the synchronization cluster analysis methods pro- 
posed in Refs. [35, 36] particularly with respect to field 
applications using model systems. In addition, we show 
that the method proposed in Ref. [36] allows to detect 
multiple synchronization clusters in long-lasting multi- 
channel EEG recordings from epilepsy patients. 

This paper is organized as follows. Since both methods 
are based on the mean phase coherence we first recall its 
definition and interpretation as a bivariate measure for 
phase synchronization (Sec. II A). In Sec. II B we briefly 
introduce the multivariate synchronization cluster meth- 
ods, namely the mean field approach (Sec. II B 1) and the 
eigenvalue decomposition approach (Sec. II B 2). Next, 
we present our simulation studies aiming at an explo- 
ration of the limitations and a comparison of both meth- 
ods (Sec. Ill A). Finally, in Sec. IIIB we present findings 
that were obtained from a spatiotemporal synchroniza- 
tion cluster analysis of multichannel EEG data recorded 
from epilepsy patients using the eigenvalue decomposi- 
tion approach. 



II. METHODS 

A. Measuring phase synchronization 

Phase synchronization was first described by Christi- 
aan Huygens [39] in the 17th century and can be defined 
as the locking of the phases of two oscillating systems j 
and k: 



A(p jk (t) = (pj(i) - ipkif) = const. 



(1) 



In a statistical way the degree of phase synchroniza- 
tion can be quantified by measuring the phase differences 
Aifjk n times and transforming them onto a unit cycle in 
the complex plane. The underlying circular distribution 
of the sample can be characterized by means of direc- 
tional statistics [40] with the mean phase coherence Rjk 
[41, 42] 
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where cpjm denotes the phase of system j in measure- 
ment m. By definition, Rjk is confined to the interval 
[0,1] where Rjk = 1 indicates fully synchronized sys- 
tems. In field applications the sample size n is typically 
limited making Rjk an estimate of the true population 
value pjk = |^e^ - ^ fc ))| of the underlying distribution 
of phase differences (the angle brackets denote the aver- 
age over all members of the population). 

When analyzing real valued time series s(t) different 
methods can be used to extract phase information. Meth- 
ods based on the Fourier, the Hilbert, or the wavelet 
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transform were shown to be equivalent under relatively 
general assumptions [43, 44]. The main idea is to map 
the data onto the complex plane using a function z and to 
take the complex argument in order to obtain the phase 
if = arg(z(£)). We here followed an approach based on 
the Hilbert transform using the analytic signal [45, 46] 



z(t) = s(t) + iHs(t) 



(3) 



where Hs(t) denotes the Hilbert transform of the signal 

s{t) 
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and V denoting the Cauchy principal value of the inte- 
gral. Application of the convolution theorem turns the 
last equation into 



Hs(t) 



T\s(t)\ • sgn{uj) 



(5) 



where T denotes the Fourier transformation and T~ x 
the inverse transformation, respectively. Thus, the imag- 
inary part of the analytic signal is obtained by shifting 
each frequency component of the original signal by it/ 2 
separately. It is important to keep in mind that phases 
may not have a physical meaning for arbitrary signals. If, 
however, the dynamics exhibit oscillations with a single 
main rhythm, then the phases are typically well defined 
[47, 48]. 

In our applications (see Sec. Ill) we computed Rjk 
from the discrete time series Sj(t) and Sk(t) as follows. 
First, the data were normalized to zero mean, which cor- 
responds to setting the DC Fourier coefficient to zero. 
Second, in order to avoid edge effects the first and last 
10% of a data window of size n' were tapered using a co- 
sine half wave (Hanning window) before performing the 
Fourier transform. Third, since the computation of the 
Hilbert transform requires integration over infinite time, 
which cannot be performed for a window of finite length, 
10% of the calculated phase time series (fj(t) and cpk(t) 
were discarded on each side of the window, reducing the 
size of the considered phase time series to n = 0.8 • n' . 
For a data window of size n and m G {0, . . . , n — 1} de- 
noting the data point within the window the mean phase 
coherence Rjk between the two time series with sampling 
interval At was obtained by identifying (pj m = cpj (mAt) 
and applying Eq. (2). 



B. Synchronization cluster analysis 

For N oscillating systems the pairwise computation of 
the mean phase coherence Rjk leads to a matrix R which 
is symmetric due to definition (2). Subsets of oscillating 
systems can be interpreted as synchronization clusters if 
these systems exhibit higher mean phase coherence val- 
ues between each other than with systems not included in 
the same subset. Allefeld and colleagues [35, 36] recently 



proposed two different multivariate approaches that aim 
at identifying synchronization clusters in which the oscil- 
lating systems participate with different strength. Both 
approaches are based on the bivariate mean phase coher- 
ence and shall briefly be recalled here. 



1. Mean field approach 

In Ref. [35] a mean field approach has been presented 
which assumes the existence of a single synchronization 
cluster C in the data. Generating a common rhythm all 
oscillating systems constitute the cluster but contribute 
to its emergence to a different extent. A mean field with 
phase <£> can then characterize the dynamics of the collec- 
tive behavior. Using Eq. (2) the degree of participation 
{participation strength) of each oscillating system j to the 
cluster C can be quantified by 



R 



jc 



1 n— 1 



Aw 



(6) 



A straightforward derivation of the phase has been 
demonstrated in Ref. [35] for a simple model system. 
This method, however, requires exact knowledge about 
the underlying equations of motion and thus can not be 
applied to unknown systems in general. Nevertheless, 
the participation strength (6) can be estimated directly 
as follows. 

If the assumption of the existence of a single cluster 
holds, a mean field can be introduced such that the dy- 
namics of the phase differences are decoupled. If, in addi- 
tion, the noise affecting the oscillating systems is statis- 
tically independent for each system the phase differences 
Aipj = (fj — <fr become independent random variables. 
Hence the population values pjk of the mean phase co- 
herences Rjk turn into 



Pjk 
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i(Aipj-A<pk) 
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Pjc ■ Pkc for j ^ k (pjj = 1). 



(7) 



Taking into account that the mean phase coherence is 
asymptotically normally distributed Rjk ~ ^{pjk-,°~jk) 
[49] and is an empirical estimate of pjk = Pjc pkc (see 
Sect. II A), a maximum likelihood estimation of the pjc 
leads to the minimization of the sum of squared weighted 
errors that defines the cost function 



j>k 



with 
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'jk — PjC PkC 
(Jjk 



(8) 



Assuming that the circular distribution of the phase dif- 
ferences can be sufficiently approximated by a wrapped 
normal distribution, the standard deviation Ojk of the 
sampling distribution of Rjk can be expressed in terms 
of p jk (cf. [40,49]) 



Ojk 



(1 



p]c Pkc)- 



(9) 



In the following we refer to pjc as defined above as Rjc, 
thereby emphasizing an interpretation as a to- cluster syn- 
chronization strength analogous to Eq. (6). In our appli- 
cations participation strengths were computed by mini- 
mizing r using an iterative algorithm proposed in Ref. 
[50]. 



2. Eigenvalue decomposition approach 

In Ref. [36] another approach has recently been pro- 
posed that is based on the eigenvalue decomposition of 
the matrix R and appears to allow identification of mul- 
tiple clusters. The procedure makes use of findings by 
Miiller and colleagues [29] who demonstrated that infor- 
mation about the correlation structure of multivariate 
data sets is imprinted into the dynamics of the eigenval- 
ues and into the structure of the corresponding eigenvec- 
tors by nonrandom level repulsion. 

The eigenvalues A c and eigenvectors V c of the symmet- 
ric and real-valued matrix R (see Eq. 2) can be obtained 
by solving the eigenvalue equation 



R 



A, 



ce{l,...,7V} 



(10) 



which, in general, has TV different solutions. In the fol- 
lowing it is assumed, that the eigenvectors are normalized 
(|z? c | = 1). Being transformed by an orthogonal transfor- 
mation into the basis of its eigenvectors R becomes the 
diagonal matrix D of its eigenvalues. The invar iance of 
the trace under this transformation leads to the equation 
N = tr(R) = tr(D) = S ^2 C A c . In the case of systems j, k 
showing no phase synchronization (Rjk = for j ^ k) the 
equation is trivially fulfilled by A c = 1 for all c, whereas 
the occurrence of entries Rjk > for j ^ k induces a level 
repulsion, a combined increase and decrease of eigenval- 
ues in such a way that TV = A c still holds. 

A multivariate analysis is realized as follows (cf. Ref. 
[36]): Each eigenvalue A c > 1 is associated with a syn- 
chronization cluster and quantifies its strength within the 
data set. The internal structure of cluster c is described 
by the corresponding eigenvector v c . Being normalized 
yC ]c — 1) ^ s components quantify the relative in- 
volvement of each system j to cluster c by i/j c . Combin- 
ing the eigenvalue A c and the index v 2 - c the "absolute" 
involvement of a system j in a cluster c can be described 
by the participation index 



Pj,c 



(11) 



Consequently, system j is considered as belonging to clus- 
ter c for which its participation index becomes maximal. 

In case of non- vanishing entries of R between systems 
belonging to different clusters (inter-cluster synchroniza- 
tion) it can be observed that these clusters are not char- 
acterized by separate but by a superposition of eigenvec- 
tors. In order to adjust the interpretation of the partic- 
ipation indices as the degree of involvement of a system 
within one cluster (see. Eq. 11) it was proposed in Ref. 



[36] to compute the in a first step and to assign the 
systems to the clusters as mentioned above. In a second 
step the matrix entries representing inter-cluster synchro- 
nization are set to zero and the participation indices are 
computed on the trimmed matrix again. In our applica- 
tions we followed exactly this scheme. 

Summarizing this section we conclude that both meth- 
ods, the mean field and the eigenvalue decomposition 
approach, seem to provide similar results in single clus- 
ter configurations for which an almost functional depen- 
dency of R 2 - c = pj y i (Ai denoting the largest eigen- 
value) can be observed [36]. Nevertheless, the question 
whether the approaches provide meaningful results in 
case of multi-cluster configurations - particularly with 
regard to field applications - remains unaddressed and 
shall be investigated in the following section. 



III. APPLICATIONS 

A. Simulations 

We here studied a three-cluster configuration which 
consisted of a lattice of 32 coupled identical Lorenz sys- 
tems (cf. Fig. 1). Each system j is defined by the differ- 
ential equations 



8 

X 3 ~ 3 X: 



+ yjZj+ej ■ (x Dl ^ 3 -Xj) (12) 



28 -z 



10- (Vj ~ z o) 



where ej denotes the coupling strength controlling the in- 
fluence of the diffusive unidirectional coupling on system 
j by one of three chosen driving systems ^1,2,3- 




FIG. 1: Three-cluster configuration consisting of a lattice of 
32 coupled Lorenz systems. Systems within a gray-shaded 
area belong to one cluster and are driven by one of the driving 
systems D\ — 10, D2 = 15 and D3 = 28 which are highlighted 
in white. 

In order to control this configuration by a single param- 
eter we set €j = e for systems within each cluster while 
€j = for the uncoupled systems. Taking randomly cho- 
sen points in the state space near the Lorenz attractor as 
initial conditions the differential equations were iterated 
using a fourth order Runge-Kutta algorithm [51] with a 
step size of 0.01. In order to eliminate transients, the first 
10 4 iterations were discarded. For increasing coupling 
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FIG. 2: (Color online) (a) Dependence of the cost function T on the coupling strength e for a single-cluster and a three- 
cluster configuration, (b) Spatial distribution of color-coded participation strength Rjc for selected values of e (see Fig. 1 for 
arrangement of clusters). 
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FIG. 3: (a) Dependence of participation strength Rjc on the 
coupling strength e for exemplary systems j G {1,8, 29}. (b) 
Dependence of cluster strength Sd,d — 1, . . . , 3 on the cou- 
pling strength e. 



strengths e = 0.0 to 1.4 (step size 0.05) we generated 
scalar time series of the x-components (n f = 5 • 10 5 data 
points) and computed the participation strength Rjc and 
the participation indices p J?c . 

For the mean field approach the configuration repre- 
sents a violation of the single-cluster assumption. When 
increasing the coupling strength e the mean phase coher- 
ence values between the coupled systems increased and 



consequently the three synchronization clusters emerged. 
This is reflected by the dependency of the cost function 
r on e as shown in Fig. 2a. Since the single-cluster as- 
sumption leading to Eq. (7) was violated, T increased 
rapidly for e > 0.4. For comparison, we repeated the 
analysis for a single-cluster configuration that was gen- 
erated by setting the coupling terms between the sys- 
tems of clusters driven by D\ and D2 to zero. Here 
T < 0.6 for all e values (cf. Fig. 2a). Using Y as an 
indicator that reflects the violation of the single-cluster 
assumption (cf. Ref. [36]), which increases with the in- 
ternal degree of synchronization between systems within 
a cluster, we generated plots (cf. Fig. 2b) showing the 
spatial distribution of the participation strength Rjc for 
the three-cluster configuration chosen here. For e < 0.4 
(r < 0.6) no clear-cut cluster structure could be identi- 
fied in the spatial distribution of participation strengths. 
For e > 0.5 (r ^> 0.6) either a single or all three syn- 
chronization clusters emerged, however, with a varying 
degree of visibility. Interestingly, for various coupling 
strengths (Fig. 2b, e.g. e e {0.6,0.9,1.1}) the partic- 
ipation strengths of systems involved in one cluster ex- 
hibited higher values than the remaining participation 
strengths. This behavior is shown in more detail in Fig. 
3a for exemplary systems involved in one cluster. In or- 
der to further elucidate this phenomenon we calculated 
the sum Sd = ^2j > k-,j,kec d Rjk, where Cd denotes the 
index set of systems belonging to cluster d. Sd quanti- 
fies the strength of cluster d and increased, on average, 
with an increasing coupling strength e (cf. Fig. 3b) while 
fluctuations can be attributed to the mean phase coher- 
ence values being asymptotically normally distributed. 
A synchronization cluster d became visible in the spa- 
tial distribution of participation strengths Rjc when its 
strength Sd dominated the configuration (cf. Fig. 3). 
This effect was caused by the cost function T exhibiting 
several local minima whose values changed according to 
the cluster strengths Sd- For e = 0.8 all Sd attained sim- 
ilar values. Here the minimization algorithm terminated 
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FIG. 4: (Color online) Spatial distribution of color-coded 
participation indices pj jC for selected values of the coupling 
strength e. Numbers denote the clusters c to which each sys- 
tem was assigned by the eigenvalue decomposition method. 
Eigenvalues were sorted in a descending order. 



because of having found a single minimum which would 
lead to the (mis) interpretation of all systems belonging 
to one single cluster. 

In contrast to the mean field method, the eigenvalue 
decomposition approach provides not just a single set of 
scalars but multiple sets of scalars representing different 
clusters. The eigenvalues were labeled here according to 
Ai > . . . > Xn thereby enabling an easy identification of 
the dominating cluster structures by the corresponding 
index c. As shown in Fig. 4 for low coupling strengths 
(e < 0.3) the method was not able to detect the three 
different clusters but indicated a multitude of small clus- 
ters. This can be attributed to random fluctuations of the 
mean phase coherences. Increasing slightly the coupling 
strength to e = 0.4 two of the three simulated clusters 
were already detected in the two largest eigenvalues Ai 
and A2 whereas the corresponding participation strengths 
did not reveal any cluster structure (cf. Fig. 2b). For 
e > 0.5 all three clusters were clearly visible and identi- 
fied by the three largest eigenvalues. Their participation 
indices increased with increasing coupling strength e (see 
e.g. e = 1.1 to 1.4) thereby reflecting the internal de- 
gree of synchronization within the clusters. The cluster 
indices c changed for different coupling strengths e due 
to the fact that the cluster strengths fluctuated (cf. Sd 
in Fig. 3b) which is reflected by the eigenvalues (see e.g. 
e G {1.1, 1.3} in Fig. 3b and Fig. 4). 

In contrast to the mean field approach the eigenvalue 
decomposition approach was capable of distinguishing 
between different clusters (see [36]). Limitations of the 
method, however, that are related to parameters influenc- 
ing the detectability of different clusters have not been 
studied so far. To address this issue we designed a con- 
figuration consisting of N = 32 systems that formed two 
clusters C\ = {1, . . . , r} and C2 = {r + 1, . . . , N} where 
the parameter r controls the relative size of an individ- 
ual cluster. We did not study dynamical systems here 
but instead assigned values to the entries of matrix R 
that were drawn from normal distributions given that 



the mean phase coherences are asymptotically normally 
distributed (see section IIB1). Mean phase coherences 
between systems of different clusters were drawn from 
the inter-cluster distribution M{p^ nt \ <j( mt )), while those 
representing the entries between systems of the same 
cluster C\ and C2 were drawn from Af(p^\ cr^) and 
M{p (y2 \(J^)^ respectively. For a given population value 
p the standard deviation a was determined by equation 
(9) using a sample size of n = 200. 

Recalling that the detection of different clusters is 
based on the participation indices computed in the first 
processing step of the method (see Sec. II B 2), which 
assigns a system to a cluster for which its participation 
index becomes maximal, we here quantify an erroneous 
assignment by considering the differences of the partic- 
ipation indices obtained from the first processing step 
X312 — Pj,i ~Pj,2 of system j and clusters C\ and C2. For 
the configuration considered here we define the weighted 
amount of erroneous assignments as 
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A 2 : N x =2 
: else 



(13) 



where N\ denotes the number of eigenvalues A > 1 and 
Ai and A2 are determined by 



A i = E e te) Ex*™ • e (^) ( 14 ) 

(N \ _1 N 

j=r+l ) j=r+l 

where is the Heaviside step function. If system j is 
erroneously assigned to the cluster C\ but belongs to a 
cluster C2 by construction, 0(xji2) wm be larger than 
zero. Thus A2 (Ai) sums up incorrect assignments of 
systems j to cluster C\ (C2) weighted by the differences 
of the corresponding participation indices \. If all sys- 
tems are assigned to the correct cluster, then A = 
by definition. Since the eigenvalues quantify the cluster 
strengths, a variation of the relative size of an individ- 
ual cluster (by varying r) makes the largest eigenvalue 
traverse a minimum. This enables us to label the eigen- 
values A c > 1, c G {1,2} according to their clusters C\ or 

c 2 . 

We studied a transition from a configuration that con- 
sisted of two clusters C\ and C2 (p^ = p^ = 0.8) to 
a single-cluster configuration Co = {1, . . . , N} by succes- 
sively increasing p( m ^ from 0.0 to 0.8 in steps of 0.01. 
Fig. 5 shows the weighted amount of erroneous assign- 
ments A depending on the relative size of an individual 
cluster r and on the inter-cluster synchronization level 



Jint) 



The eigenvalue decomposition method success- 



fully identified (A = 0) the two clusters for low values 
of the inter-cluster synchronization level (p( m£ ) < 0.17) 
except for clusters of equal size (r = N/2). The weighted 
amount of erroneous assignments rapidly increased, how- 
ever, for increasing p^ int \ thereby failing to detect the 
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FIG. 5: Dependence of the error A on the relative size r of 
an individual cluster and on the inter-cluster synchronization 
level p (mt) . 



(b) 




FIG. 6: Dependence of the error A (a) and N\ (b) on the 
relative cluster size r and on the mean cluster synchronization 
level // 2 ) of cluster Ci. 



cluster structures for an increasing set of r values. The 
number of eigenvalues N\ being larger than 1 appeared to 
be a sensitive measure for the number of clusters present 
in the data. We observed N\ — ^ 1 (i.e., A = by defi- 
nition, see Eq. 13) only for relatively large values of the 
inter-cluster synchronization level (p( mt ) > 0.7) where 
the transition to a single-cluster configuration was com- 
pleted. 



Next we studied a transition from a configuration of 
two clusters d and C 2 = p (2) = 0.8, p^ = 0.2) 
to a configuration with only one cluster (here C\) by de- 
creasing // 2 ) in steps of 0.01 down to 0.2 (see Fig. 6a). 
The weighted amount of erroneous assignments A van- 
ished for // 2 ) approaching the mean synchronization level 
of the inter-cluster distribution (p( mt ) = 0.2). In this 
range the number of eigenvalues N\ increased (Fig. 6b). 
These additional structures were formed by those sys- 
tems, which were previously involved in cluster C2 for 
higher values of p^ 2 \ This fragmentation of a cluster c 
can also be observed when increasing the corresponding 
standard deviation independently from the popula- 
tion value p^ c \ thereby taking into account various un- 
certainties (e.g. definition of observables, measurement 
precision, finite sample size, phase extraction). Since the 
additional structures were caused by random fluctuations 
of the mean phase coherences within the cluster, we re- 
fer to them as pseudo- clusters in the following. In our 
simulations we observed pseudo-clusters of different sizes, 
whose corresponding eigenvalues were only slightly larger 
than unity compared to the eigenvalues of the true clus- 
ters. 

When varying p^ the set of r values for which A > 
remained almost stable but was shifted to lower r val- 
ues (Fig. 6a). Therefore, together with the relative size 
of an individual cluster the level of synchronization be- 
tween systems within a cluster must be regarded as cru- 
cial for the eigenvalue decomposition approach. The re- 
sults shown in Figs. 5 and 6a can then be sufficiently 
explained with the help of the cluster strength Sd which 
considers both, the cluster size and the level of synchro- 
nization between systems within a cluster. The method 
failed to correctly assign all systems to their clusters in 
case of similar cluster strengths \S± — S2I < ft, where 
n increases with an increasing level of inter-cluster syn- 
chronization. The latter leads to superpositions in the 
eigenvector components and, as a result, makes it im- 
possible for the method to correctly detect the clusters 
by the maximum participation index criterion in certain 
configurations. In order to further elucidate this phe- 
nomenon we show in Fig. 7 the components Vj C of the 
eigenvectors corresponding to the two largest eigenvalues 
for the case r = 16 and p^ = 0.8. The squared compo- 
nents i/j C of both vectors attained similar values, causing 
the method to classify all systems as belonging to a single 
cluster. Nevertheless, the two clusters C\ and C2 were 
clearly visible in the linear combinations V\ ±1*2. 

Summarizing this section we conclude that the eigen- 
value decomposition approach successfully identified the 
three synchronization clusters in our simulations us- 
ing coupled Lorenz systems and, moreover, appears to 
be sensitive to detect clusters even for low coupling 
strengths. Although the mean field approach is by defi- 
nition not able to distinguish between different clusters, 
the method appears to be relatively robust even in situa- 
tions where the single-cluster assumption is not (entirely) 
fulfilled, leading to higher participation strengths for sys- 
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FIG. 7: Components j of eigenvectors belonging to the two 
largest eigenvalues Ai and A 2 for the case of r = 16, = 
p (2) = 0.8, p (tnt) = 0.2. 



terns of the dominating cluster in our simulation. When 
comparing both methods the eigenvalue decomposition 
approach can be regarded as superior when analyzing 
data exhibiting different synchronization clusters, which 
can be expected e.g. for multichannel EEG data. 



B. EEG data 

In this section we present exemplary findings obtained 
from applying the eigenvalue decomposition method to 
multichannel EEG time series that were recorded from 
three epilepsy patients (denoted as PI, P2, and P3) suf- 
fering from pharmacoresistant focal epilepsies of neocor- 
tical origin. For these patients complete seizure control 
can be obtained by resecting the part of the brain re- 
sponsible for seizure generation (epileptic focus). This 
requires an exact localization of the epileptic focus and 
its delineation from functionally relevant brain structures 
during the presurgical workup. When no concordant in- 
formation can be achieved from noninvasive diagnostic 
techniques, the EEG is recorded from implanted elec- 
trodes over a longer period, typically 2-3 weeks. The 
analyses reported here were made after surgery had taken 
place, and after it had become clear from its success 
whether the localization of the epileptic focus had been 
correctly predicted. All patients had signed informed 
consent that their clinical data might be used and pub- 
lished for research purposes, and the study was approved 
by the local medical ethics committee. 

Previous studies have shown that even during seizure- 
free intervals the seizure generating area of the brain ex- 
hibited higher interdependences [52] and an higher degree 
of synchronization [42] than other brain areas. Together 
with results obtained from applying univariate time se- 
ries analysis techniques (see e.g. [53] and references 
therein) these findings allow an improved understand- 
ing of intermittent dysfunctioning of the brain between 




FIG. 8: Schematic view of the electrode grid with N = S x 4 
contacts placed over the left temporal lateral neocortex. 



seizures and provide potentially useful diagnostic infor- 
mation. We here addressed the question whether com- 
plementary information can also be obtained from a mul- 
tivariate approach. Specifically, we investigated whether 
the eigenvalue decomposition approach (a) allows to lo- 
calize the epileptic focus analyzing EEG recordings from 
the seizure-free interval only, and (b) is capable of detect- 
ing short time changes of synchronization patterns asso- 
ciated with physiological processes in the human brain. 
To this end we analyzed continuous EEG recordings that 
lasted 39 h for patient PI and 26 h for patient P2 cover- 
ing different physiological and pathophysiological states 
of the patients. In addition, we analyzed an EEG record- 
ing lasting 30 minutes during which patient P3 was si- 
multaneously presented item pairs, either word pairs or 
unpronounceable letter string pairs of 5-11 letters length, 
and was instructed to press a button if he/she recognized 
word pairs with identical or highly similar meanings (syn- 
onyms) (task Tl) or identical letter strings (task T2). 
No button press was demanded in case of non-synonym 
word pairs or letter strings in which one consonant dif- 
fered. Tasks alternated every 3 min. The whole sequence 
was adopted from Ref. [54] and comprised three blocks of 
word and letter string pairs each. A baseline recording of 
10 minutes was performed with eyes open before and af- 
ter the experiment. Recent findings [55, 56] indicate that 
tasks such as T2 are associated with a greater demand for 
decision-making processes due to the involvement of dif- 
ferent phenomena like reading, phonological retrieval, or 
orthographic analysis during such conditions. We thus 
hypothesized to observe a higher degree of collectivity 
among neuronal assemblies particularly in brain areas as- 
sociated with language processing. We did not expect to 
observe synchronization phenomena associated with the 
epileptogenic process in these brain structures from pa- 
tient P3 since the epileptic focus was localized in a more 
distant brain structure. 

The EEG was measured from grid electrodes (rectan- 
gular flexible grids of N = 8 x 4 contacts) placed onto 
the temporal lateral neocortex (see Fig. 8). EEG data 
were sampled at 200 Hz using a 16-bit analog-to-digital 
converter and filtered within a frequency band of 0.3-70 
Hz. Using a moving- window technique EEG signals were 
divided into segments of n r = 4096 sampling points each, 
and segments overlapped by 20%. The length of the re- 
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time 30 min 



FIG. 9: Examples of the spatio-temporal evolution of the 
cluster defined by the largest eigenvalue. EEG data from pa- 
tient PI recorded during day-time (top) and during night-time 
(bottom). For each segment (duration: 16.38 s) channels be- 
longing to the cluster are drawn in black. Note that channels 
Al and A2 were used as reference during the recording. 



suiting segments corresponded to 16.38 s at the given 
sampling rate and can be regarded as a compromise be- 
tween the required statistical accuracy for the calculation 
of the mean phase coherence and the approximate sta- 
tionarity within a segments length (see Ref. [57] and 
references therein). After the calculation of the mean 
phase coherence (cf. Sect. II A) for all channel combina- 
tions from each segment, the eigenvalue decomposition 
method was applied (cf. Sect. II B 2). 

In order to estimate the number of clusters being de- 
tected on average in the data we sorted eigenvalues and 
corresponding eigenvectors of each segment in a descend- 
ing order (Ai > . . . > Ajv) and computed the mean val- 
ues A c over all segments for each patient. The five to 
six largest averaged eigenvalues A c were larger than 1 
causing the method to detect on average an according 
number of clusters per segment. Based on our findings 
presented in Sect. Ill A we expected the lower part of 
the eigenvalue spectrum being larger than 1 to yield a 
considerable amount of pseudo-clusters. We therefore re- 
stricted further analyses to the three largest eigenvalues. 

Since a cluster is represented by a set of participating 
channels, each cluster can be written as a bit string of size 
N where the bits represent the channels being involved 
(1) or not involved (0) in the cluster. This allowed us to 
discard the information about the relative involvement of 
each system within one cluster as reflected by the partic- 
ipation indices and to handle the clusters in a convenient 
way, namely by considering their simplified representa- 
tion. Bit strings connected to eigenvalues fulfilling the 
threshold criterion (A > 1) varied largely over time as 
shown exemplar ily in Fig. 9. Apart from frequently ap- 
pearing bit strings in which all channels were involved, 
structures in time could be observed, where successive 
bit strings differed only in few bits. Given that the mean 
phase coherence fluctuates (due to measurement noise, a 
limited number of data points, or even due to short-term 
physiological or pathophysiological phenomena within a 
segment), the same cluster structure cannot be expected 
to show up in exactly the same bit string representation 
but may vary in some bits. In order to minimize this ef- 
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FIG. 10: Drawings of the grid with electrode contacts Al (top 
left), A8 (top right) to D8 (bottom right). Channels partic- 
ipating in a cluster are marked by the same letter. Because 
of the averaging applied here a channel can be observed to 
participate in multiple clusters for different segments. 



feet we sorted bit strings that were visible for more than 
6 minutes in the data into groups which differed only 
in up to three bits. Moreover, we discarded groups of 
bit strings which represented clusters of size 0, 1, or 32 
channels. In Fig. 10 we show groups that represented 
the largest number of bit strings in the EEG data. Inter- 
estingly, for both patients the spatial distribution of elec- 
trodes involved in the cluster labeled "a" corresponded 
to the spatial extent of the epileptic focus as determined 
by the presurgical workup. This cluster could be ob- 
served throughout the datasets and in total for at least 
112 minutes for patient PI and 125 minutes for patient 
P2. Indeed, the patients were operated on exactly this 
region and are now free of seizures. Thus, we conclude 
that the method seems to be quite sensitive for detect- 
ing synchronization clusters in EEG time series recorded 
from epilepsy patients even during the seizure-free inter- 
val. 

Whereas the aforementioned results were obtained by 
analyzing the bit string groups which corresponds to a 
temporal average (over all segments), we now investi- 
gated whether the eigenvalue decomposition approach al- 
lows to detect short-term changes of cluster structures 
that can be related to physiological synchronization phe- 
nomena (the language processing paradigm mentioned 
above). The participation indices were computed from 
the EEG recording from patient P3 and translated into 
bit strings using the same criteria as mentioned above. 
Bit strings corresponding to the three largest eigenval- 
ues are shown in Fig. 11. Interestingly, a modulation of 
the bit strings depending on the tasks can be observed. 
The occurrence of bit strings representing a cluster of 
the electrode contacts A7, B7-8, and C7-8 is noticeable 
particularly during task T2 (letter matching task). For 
this patient the brain area covered by these electrode 
contacts was associated with language processing (Wer- 
nicke's area). The observed higher level of synchroniza- 
tion probably reflects the higher degree of collectivity 
among neuronal assemblies that are involved in the cog- 
nitive operations comprising this task. This suggests that 
the eigenvalue decomposition method is capable of de- 
tecting short time changes of synchronization patterns as- 
sociated with physiological processes in the human brain. 
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FIG. 11: (Color online) Spatio-temporal evolution of clusters denned by the largest, second-, and third- largest eigenvalues 
calculated from EEG data from patient P3 recorded during a language processing paradigm involving two different tasks 
(separated by black lines). The patient was simultaneously presented item pairs, either word pairs or unpronounceable letter 
string pairs, and was instructed to press a button upon recognizing synonyms (Tl) or identical letter strings (T2). A baseline 
recording (B) was performed before and after the experiment. Each pixel in the figure represents the color-coded cluster 
membership of a given electrode contact (A1-A8, B1-B8, C1-C8, D1-D8 shown on the ordinate; cf. Fig. 8) and for a given 
EEG segment of duration 16.38 s (abscissa). Contacts A7, B7-8, and C7-8 covered a brain area associated with language 
processing (Wernicke's area) and were most noticeable during task T2. 



IV. CONCLUSION 

In this paper we have studied two multivariate phase 
synchronization analysis methods, namely a mean field 
approach [35] and an eigenvalue decomposition approach 
[36]. While the mean field approach assumes the exis- 
tence of a single synchronization cluster in the data, the 
eigenvalue decomposition approach appears to be capable 
of identifying multiple clusters. Based on the results of 
numerical simulations of multiple synchronization clus- 
ters within a lattice of 32 coupled identical Lorenz sys- 
tems, we demonstrated that the mean field approach ap- 
pears to be relatively robust even in situations where the 
single-cluster assumption is not entirely fulfilled. The 
eigenvalue decomposition approach successfully identi- 
fied multiple synchronization clusters in our simulations 
and appears to be sensitive to detect clusters even for 
weak couplings. However, in case of non- vanishing inter- 
cluster synchronization the method failed to correctly 
assign the systems to their clusters in certain configu- 
rations. The influence of measurement noise, which was 
not discussed in depth here, can be expected to lead to 
the fragmentation of clusters present in the data. Nev- 
ertheless, the eigenvalue decomposition approach can be 
regarded as superior when analyzing data exhibiting dif- 
ferent synchronization clusters. 

When being applied to field data a number of influenc- 



ing factors limit the significance of the eigenvalue decom- 
position approach and need further investigations. These 
factors include the finite length of available data, random 
spatio-temporal correlations, or even non-random corre- 
lations being induced by the data acquisition system (e.g. 
filtering or, as in the case of EEG recordings, the choice of 
a suitable reference). Nevertheless, our preliminary ap- 
plications of the eigenvalue decomposition approach to 
multichannel EEG recordings from epilepsy patients in- 
dicate that the method allows to gain deeper insights into 
the collective dynamics of neuronal networks, both under 
physiological and pathophysiological conditions. Despite 
the limited EEG database used in this study, the achieved 
results can be regarded as promising. Further evaluations 
on a larger EEG database are currently underway. 
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